#include "Pythia.h"
#include "TH1D.h"
#include "TH1D.h"
#include <vector>
#include <iostream>
#include "TROOT.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TRandom.h"
#include "TF1.h"
#include "TFile.h"
#include "TChain.h"
#include "TClonesArray.h"
#include "TLorentzVector.h"
#include <TStopwatch.h>
#include "Case_Table.h" // all const valuse and setting is here



using namespace std; // 
using namespace Pythia8; //  

typedef vector<double> DOUBLE1D;
typedef vector<DOUBLE1D> DOUBLE2D; // define new variable type for 2Dim vector
typedef vector<DOUBLE2D> DOUBLE3D; // define new variable type for 3Dim vector

int main(int argc, char **argv)
{

  TROOT root("flow_jet","run mc");
  if ( argc<6 ) {
    cout<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
    cout<<"usage :  "<<argv[0]<<" <outputFile> <N_ch TPC> <N_ch V0A> <N_ch V0C> <Rndm_seed> <Nevt>"<<endl;
    cout<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
    cout << endl << endl;
    exit(1);
  }

  // CONSTANT
  Int_t nPar[ND];
  char *outFile = argv[1];
  nPar[kTPC]  = atoi(argv[2]); // vary with multiplicity
  nPar[kV0A]  = atoi(argv[3]); // vary with multiplicity
  nPar[kV0C]  = atoi(argv[4]); // vary with multiplicity
  int Rndm_seed = atoi(argv[5]);

  if(ND>3){
    nPar[kjTPC] = nPar[kTPC];
    nPar[kjV0A] = nPar[kV0A];
    nPar[kjV0C] = nPar[kV0C];
  };

  Int_t numberEvents= atoi(argv[6]);
  //char *cardName = argv[8]; // input card file name

  TFile *fout = new TFile(outFile,"recreate");

  cout<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
  cout << "Outfile : "<< outFile << endl;
  cout << "Inputs : "<<endl;
  cout <<"\t nParTPC = "<< nPar[kTPC] << endl;
  cout <<"\t nParV0A = "<< nPar[kV0A] <<endl;
  cout <<"\t nParV0C = "<< nPar[kV0C] <<endl;
  cout<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;

  // ---------- prepareing histograms --------------
  // ******************** Histograms for Detector information ********************************
  //for detector
  TH1D *hphi[ND];
  TH1D *hphi_m[ND];
  TH1D *hDeltaphi_measure[ND];
  TH1D *hJet_phi[ND];
  TH1D *hNjetparticles[ND]; 
  for(int id=0; id<ND; id++){
    hphi[id] = new TH1D(Form("hPhi_%d",id),Form("Phi_%d",id), 100, -pi, pi);
    hphi_m[id] = new TH1D(Form("hPhi_m_%d",id),Form("Measured Phi_%d",id), 100, -pi, pi);
    hDeltaphi_measure[id] = new
        TH1D(Form("hDeltaphi_measure_%d",id),Form("hDeltaphi_measure_%d",id),100,-2*pi,2*pi);
    hJet_phi[id] = new
        TH1D(Form("hJet_phi_%d",id),Form("phi distribution of jet particles in %s"
                                         ,vdets_str[id].Data()),100,-2*pi, 2*pi);
    hNjetparticles[id] = new
        TH1D(Form("hNjetparticles_%d",id),Form(" Njet particles in %s",
                                               vdets_str[id].Data()),100,0,100);
  };

  // *****************************************************************************************
  // ******************** Histograms for Event plane information ********************************
  TH1D *hPsi[ND][NH];
  TH1D *hPsi_m[ND][NH];
  TH1D *hdpsi[ND][NH];
  TH1D *hCosndpsi[ND][NH]; // for calculate resolution resolution
  TH1D *h2sub_Cosndpsi[ND][NH]; // for calculate resolution resolution with 2sub
//  TH1D *hEPResolution[ND][NH];
  TH1D *hJetEPdPhi[ND][NH];
  TH1D *hJetEPCosdPhi[ND][NH];
  TH1D *hEP_pTLphi[ND][NH];
//  TH1D *hEPResolution_2sub[ND][NH];
//  TH1D *hEPResolution_3sub[ND][NH];


  for(int ih=2;ih<NH;ih++) {
    for(int id=0;id<ND;id++) {
      hPsi[id][ih] = new TH1D(Form("hPsi_%02d%02d",id,ih),
                              Form("Psi distribution"), 100, -pi, pi);
      hPsi_m[id][ih] = new TH1D(Form("hPsi_m_%02d%02d",id,ih),
                              Form("Measured Psi distribution"), 100, -pi, pi);
      hdpsi[id][ih] = new TH1D(Form("hdpsi_%02d%02d",id,ih),
                               Form("Detector = %s, %d th Generated psi - Calculated psi"
                                    ,vdets_str[id].Data(),ih), 100, -pi/ih, pi/ih);
      hCosndpsi[id][ih] = new
          TH1D(Form("hCosndpsi_%02d%02d",id,ih),Form("Detector = %d, %dth Harmonics resolution",id,ih) , 1000, -1.1, 1.1);

      h2sub_Cosndpsi[id][ih] = new
          TH1D(Form("h2sub_Cosndpsi_%02d%02d",id,ih),Form("Detector = %d, %dth Harmonics resolution",id,ih) , 1000, -1.1, 1.1);

 //     hEPResolution[id][ih] = new
 //         TH1D(Form("hEPResolution_%02d%02d",id,ih),Form("Detector = %d, %dth Harmonics resolution",id,ih) , 1000, -1.1, 1.1);

      hJetEPdPhi[id][ih] = new
          TH1D(Form("hJetEPdPhi_%02d%02d",id,ih),Form("hJetEPdPhi_%s  %dth",
                                                      vdets_str[id].Data(),ih),100,-2*pi,2*pi);
      hJetEPCosdPhi[id][ih] = new
          TH1D(Form("hJetEPCosdPhi_%02d%02d",id,ih),Form("hJetEPCosdPhi_%s  %dth",
                                                         vdets_str[id].Data(),ih),100,-1.1,1.1);
      hEP_pTLphi[id][ih] = new
          TH1D(Form("hEP_pTLphi_%02d%02d",id,ih),
               Form("%s EP psi - pTL phi angle, %dth",vdets_str[id].Data(),ih),
               36, -1*pi, pi); 

//     hEPResolution_2sub[id][ih] = new
//          TH1D(Form("hEPResolution_2sub_%02d%02d",id,ih),Form("Detector = %d, %dth Harmonics ressolution",id,ih), 1000, -1.1, 1.1);
//      hEPResolution_3sub[id][ih] = new
//          TH1D(Form("hEPResolution_3sub_%02d%02d",id,ih),Form("Detector = %d, %dth Harmonics ressolution",id,ih), 1000, -1.1, 1.1);



    }
  }
  // *********************************************************************************************
  // ******************** Histograms for flow  information ****************************************
  TH1D *hEPCosndPhi[NType][NH]; // mean of this value is vn = <cos ih * dphi> : dphi= phi- psi_reco
  TH1D *hRealEPCosndPhi[NType][NH]; // vn from real EP value = <cos ih *dphi> : dphi= phi- psi_in
  TH1D *hEPVn_Evt[NType][NH];
  TH1D *hEPVn_Trk[NType][NH];
  TH1D *hRealEPVn_Evt[NType][NH];
  TH1D *hRealEPVn_Trk[NType][NH];
  TH1D *hEPdPhi_JetOnly[NType][NH];
  TH1D *hEPdPhi_BulkOnly[NType][NH];

  for(int ih=2; ih<NH; ih++){
    for(int it=0; it<NType; it++){
      hEPCosndPhi[it][ih] = new TH1D(Form("hEPCosndPhi_%02d%02d",it,ih), Form("cos%d#Delta#phi Case= %d",ih,it), 1000, -1.1, 1.1);
      hRealEPCosndPhi[it][ih] = new TH1D(Form("hRealEPCosndPhi_%02d%02d",it,ih), Form("cos%d#Delta#phi Case= %d",ih,it), 1000, -1.1, 1.1);
      hEPVn_Evt[it][ih] = new TH1D(Form("hEPVn_Evt_%02d%02d",it,ih), Form("cos%d#Delta#phi Case= %d",ih,it), 1000, -1.1, 1.1);
      hEPVn_Trk[it][ih] = new TH1D(Form("hEPVn_Trk_%02d%02d",it,ih), Form("cos%d#Delta#phi Case= %d",ih,it), 1000, -1.1, 1.1);
      hRealEPVn_Evt[it][ih] = new TH1D(Form("hRealEPVn_Evt_%02d%02d",it,ih), Form("cos%d#Delta#phi Case= %d",ih,it), 1000, -1.1, 1.1);
      hRealEPVn_Trk[it][ih] = new TH1D(Form("hRealEPVn_Trk_%02d%02d",it,ih), Form("cos%d#Delta#phi Case= %d",ih,it), 1000, -1.1, 1.1);
      hEPdPhi_JetOnly[it][ih] =  new TH1D(Form("hEPdPhi_JetOnly_%02d%02d",it,ih),
                                          Form("cos%d#Delta#phi Case= %d of Jet particles",ih,it), 400, -pi, pi);
      hEPdPhi_BulkOnly[it][ih] = new TH1D(Form("hEPdPhi_BulkOnly_%02d%02d",it,ih), 
                                          Form("cos%d#Delta#phi Case= %d of Bulk particles",ih,it),  400, -pi, pi);
    }
  }

  // *********************************************************************************************
  //******************for injected jet particle distribution ***************************************//  
  const int NBINS=200;
  double LogBinsX[201], LimL=0.1, LimH=110;
  double logBW = (log(LimH)-log(LimL))/NBINS;
  for(int ij=0;ij<=NBINS;ij++) LogBinsX[ij]=LimL*exp(ij*logBW);

  TH1D *hChargedPt= new TH1D("hChargedPt"," ", NBINS, LogBinsX ); hChargedPt->Sumw2();
  TH1D *hChargedPtDetector[ND];
  for(int id=0; id<ND; id++){
    hChargedPtDetector[id] = new TH1D(Form("hCharged%02d",id),"",NBINS, LogBinsX);
    hChargedPtDetector[id]->Sumw2();
  }
  TH1D *hChargedEta = new TH1D("hChargedEta"," ", 200, -10,10 ); hChargedEta->Sumw2();
  TH1D *hPartonPt= new TH1D("hPartonPt"," ", NBINS, LogBinsX ); hPartonPt->Sumw2();
  TH1D *hLeadingpPt= new TH1D("hLeadingpPt","",NBINS, LogBinsX); hLeadingpPt->Sumw2();
  TH1D *hLeadingpEta= new TH1D("hLeadingpEta","",200, -10, 10);  hLeadingpEta->Sumw2();
  TH1D *hLeadingpPhi= new TH1D("hLeadingpPhi","",200, -2*pi, 2*pi);  hLeadingpPhi->Sumw2();
  LimL=0.1;  LimH=300;
  logBW = (log(LimH)-log(LimL))/NBINS;
  for(int ij=0;ij<=NBINS;ij++) LogBinsX[ij]=LimL*exp(ij*logBW);
  TH1D *hJetPt    = new TH1D("hJetPt"," ", NBINS, LogBinsX ); hJetPt->Sumw2();
  //***********************************************************************************************//


  /// PYTHIA - load pythia dataset
  TChain *pythia_data = new TChain("pythiaTree");
  for( int ipythiaFile=0; ipythiaFile < 20; ipythiaFile++) {
    pythia_data->Add(Form("./pythia_data/pythia_%d.root",ipythiaFile));
  } 
  TClonesArray *pythiaTracks = new TClonesArray("TLorentzVector",1000);
  pythia_data->SetBranchAddress("tracks",&pythiaTracks);
  TLorentzVector *p = NULL;
  
  if (pythia_data->GetEntries() < numberEvents){
    cout << "Number of events are too large" << endl;
    exit(1);
  }

  int ieout = numberEvents/10;
  if (ieout<1) ieout=1;



  TRandom *myRandom = new TRandom(Rndm_seed);

  TStopwatch timer;
  timer.Start();
  //*****************************************************************************************//y
  // ---------------------  event loop start --------------- //
  //*****************************************************************************************//y

  for(int ievt =0; ievt<numberEvents; ievt++){
    if( ievt%ieout==0) {
      cout << " Event: " << ievt << " " << ievt*100.0/numberEvents <<"% "<< endl;
    }
    // print out current evnet number every 10% 


    // PYTHIA - load event from dataset  
    pythiaTracks->Clear(); //clear array first
    pythia_data->GetEntry(ievt); // load ievt in pythia dataset


     
    float pTL=0;
    int pTL_index=-1;
    double pTL_phi=-99;
  
    for( int ip=0;ip<pythiaTracks->GetEntriesFast();ip++ ){
      p = (TLorentzVector*) pythiaTracks->At(ip);
      if(p->Pt() >= pTL){ pTL= p->Pt(); pTL_index=ip; pTL_phi=p->Phi();} 
    }; //find leading particle

    p = (TLorentzVector*) pythiaTracks->At(pTL_index); // load leading particle
	double pTL_eta = p->Eta();
    //if ( pTL_eta > pTL_eta_max || pTL_eta<pTL_eta_min) continue;//cut condition
    hLeadingpPt->Fill(p->Pt());
    hLeadingpEta->Fill(p->Eta());
    hLeadingpPhi->Fill(p->Phi());

    // random event plane v0, v1 = 0
    double psi_in[NH] = {0, 0, myRandom->Uniform(-1*pi/2,pi/2), myRandom->Uniform(-1*pi/3,pi/3)};

    // v_n formular //
    TF1 *dNdphi_bulk = new TF1("dNdphi_bulk",
                               "[0]*(1+2*[1]*cos(x-[4])+2*[2]*cos(2*(x-[5]))+2*[3]*cos(3*(x-[6])))",-1*pi,pi);
    dNdphi_bulk->SetParameter(0, input_vn[0]); // underyn const(?)
    dNdphi_bulk->SetParameter(1, input_vn[1]);  //v1  - ignore directed flow
    dNdphi_bulk->SetParameter(2, input_vn[2]);  // v2 - from header file 
    dNdphi_bulk->SetParameter(3, input_vn[3]);  // v3 - from header file
    dNdphi_bulk->SetParameter(4, psi_in[1]); // v1 event plane  - 0
    dNdphi_bulk->SetParameter(5, psi_in[2]); // v2 event plane  - random generated 
    dNdphi_bulk->SetParameter(6, psi_in[3]); // v3 event palne  - random generated



    // ************************  inject bulk particle **************************
    //generatation particle with random angle phi(which followed dNdphi folmuar)
    vector <double> phis[ND];
    vector <double> phis_measured[ND]; // phis after apply detector resolution
    vector <int> phis_measured_sector[ND]; // phis sector
    for (int id=0; id<ND; id++){
      if (id == kV0AC) continue;
      if (id == kjV0AC) continue;
      for (int ip=0; ip<nPar[id]; ip++){
        phis[id].push_back(dNdphi_bulk->GetRandom());
      }
    }
    // inject bulk particles done.

    // ************************  inject jet particle **************************
    //input jet from PYTHIA
    for ( int ip=0;ip<pythiaTracks->GetEntriesFast();ip++ ){  // PYTHIA particle loop
      p = (TLorentzVector*)pythiaTracks->At(ip);

      for (int id=0; id<ND; id++){
        if (p->Eta() > Detector_eta_min[id] && p->Eta() < Detector_eta_max[id] &&
            p->Pt() > Detector_Jet_pT_min[id] && Inject_Jet[id] > 0){
          phis[id].push_back( p->Phi() );
          hChargedPtDetector[id]-> Fill(p->Pt());
          if(id == kjV0A || id== kjV0C) hChargedPtDetector[kjV0AC]->Fill(p->Pt());
          hJet_phi[id]-> Fill(p->Phi());
        }
      }
      hChargedPt->Fill(p->Pt());
      hChargedEta->Fill(p->Eta());
    }
    // calculate and save N jet particles in each detectors.
    for(int id=0; id<ND; id++){
      hNjetparticles[id]->Fill(phis[id].size() - nPar[id]); 
    }
    // injection jet done

    // *************************************************************************
    // apply detector resolution - smaering
    for (int id=0; id<ND; id++){
      for (UInt_t ip=0; ip<phis[id].size(); ip++){
        if(Detector_Resolution[id]==0)
          phis_measured[id].push_back(phis[id][ip]);
        if(Detector_Resolution[id]>0){
          //make -pi ~ pi to 0 ~ 2pi
          phis[id][ip] = phis[id][ip] + pi;
          phis_measured_sector[id].push_back(int(phis[id][ip]/Detector_Resolution[id]));
          phis_measured[id].push_back(int(phis[id][ip]/Detector_Resolution[id]) *
                                      Detector_Resolution[id] +
                                      Detector_Resolution[id]/2. );
          //make 0 ~ 2pi to -pi ~ pi
          phis_measured[id][ip] = phis_measured[id][ip] - pi;
          phis[id][ip] = phis[id][ip] - pi; 
        }
        hphi[id]->Fill(phis[id][ip]);
        hphi_m[id]->Fill(phis_measured[id][ip]);
        hDeltaphi_measure[id]->Fill( phis_measured[id][ip]-phis[id][ip] );        
      } 
    }
    // smaring done
    // *************************************************************************

    // V0AC(V0A+V0C)
    for (UInt_t ip=0; ip<phis_measured[kV0A].size(); ip++){ phis_measured[kV0AC].push_back(phis_measured[kV0A][ip]);
											phis_measured_sector[kV0AC].push_back(phis_measured_sector[kV0A][ip]); }
    for (UInt_t ip=0; ip<phis_measured[kV0C].size(); ip++){ phis_measured[kV0AC].push_back(phis_measured[kV0C][ip]); 
											phis_measured_sector[kV0AC].push_back(phis_measured_sector[kV0C][ip]); }
    for (UInt_t ip=0; ip<phis_measured[kjV0A].size(); ip++){ phis_measured[kjV0AC].push_back(phis_measured[kjV0A][ip]); 
											phis_measured_sector[kjV0AC].push_back(phis_measured_sector[kjV0A][ip]); }
    for (UInt_t ip=0; ip<phis_measured[kjV0C].size(); ip++){ phis_measured[kjV0AC].push_back(phis_measured[kjV0C][ip]); 
											phis_measured_sector[kjV0AC].push_back(phis_measured_sector[kjV0C][ip]); }




    //************************calcuate Psi with above particles************************
    // event plane calculation
    DOUBLE2D Qx( ND, DOUBLE1D( NH, 0 ) );
    DOUBLE2D Qy( ND, DOUBLE1D( NH, 0 ) );
    DOUBLE2D psi_reco( ND, DOUBLE1D( NH, 0 ) );

    for( int id = 0;id< ND; id++ ){
      for( int ih=2;ih<NH;ih++){
        Qx[id][ih] = 0;
        Qy[id][ih] = 0; // initializing

        for( UInt_t ip=0;ip<phis_measured[id].size();ip++ ){
          Qx[id][ih] += cos(ih*phis_measured[id][ip]);
          Qy[id][ih] += sin(ih*phis_measured[id][ip]);
        }
        psi_reco[id][ih] = TMath::ATan2(Qy[id][ih], Qx[id][ih])/ih; // EP_cal
        hPsi[id][ih]->Fill(psi_reco[id][ih]);

        double dpsi = psi_in[ih]-psi_reco[id][ih];

        // psi shoud beinside -pi/ih ~ pi/ih
        if(dpsi < -1*pi/ih) dpsi = dpsi + 2*pi/ih;
        if(dpsi > pi / ih ) dpsi =dpsi  - 2*pi/ih;    // trick to make hist  -pi/ih~pi/ih
        hdpsi[id][ih]->Fill( dpsi ); 
        hCosndpsi[id][ih]->Fill(cos(ih*(psi_in[ih]-psi_reco[id][ih])));   //resolution calculation

        for(UInt_t ip=nPar[id]; ip<phis[id].size(); ip++){ //loop for jet particle
          hJetEPdPhi[id][ih]->Fill( phis_measured[id][ip] - psi_reco[id][ih] );
          hJetEPCosdPhi[id][ih]->Fill(cos( phis_measured[id][ip]-psi_reco[id][ih]) );
        }
      }
    } // event plane calculration done.


    //2sub event resolution***********************************************
    double Qx_2sub[2];
    double Qy_2sub[2];
    double psi_reco_2sub[2];

    for( int id=0 ; id<ND; id++ ) { //for each detector
      if (Detector_Resolution[id] == 0 ){
        for (int ih=2 ; ih<NH; ih++ ) {  //for each harmonics
          Qx_2sub[0] = 0; Qy_2sub[0] = 0; // initializing 1st group of 2sub event
          Qx_2sub[1] = 0; Qy_2sub[1] = 0; // initializaing 2nd group of 2sub event
          for ( UInt_t ip=0 ; ip<phis_measured[id].size(); ip++){ //for each  particles
            Qx_2sub[ip%2] += cos(ih*phis_measured[id][ip]);
            Qy_2sub[ip%2] += sin(ih*phis_measured[id][ip]);
          }//end particles
          psi_reco_2sub[0] = TMath::ATan2(Qy_2sub[0], Qx_2sub[0] ) / ih;
          psi_reco_2sub[1] = TMath::ATan2(Qy_2sub[1], Qx_2sub[1] ) / ih;
          double dCos_2sub_psi = cos(ih*(psi_reco_2sub[0]-psi_reco_2sub[1]));
          h2sub_Cosndpsi[id][ih]->Fill(dCos_2sub_psi);
        } //end harmonics
      }// 2sub event for TPC done

      if (Detector_Resolution[id] != 0) {
        for (int ih=2 ; ih< NH; ih++){
          Qx_2sub[0]=0; Qy_2sub[0]=0;
          Qx_2sub[1]=0; Qy_2sub[1]=0; // initializing

          int Sub_index[V0_N_sectors]; // N_sectors are 16 or 32? 
          for (int i_index=0; i_index < V0_N_sectors; i_index++){
//            Sub_index[i_index] = myRandom->Integer(2); // get 0 or 1 for each sectors
            Sub_index[i_index] = i_index % 2; // 0, 1, 0, 1....
          }

          for (UInt_t ip=0 ; ip<phis_measured[id].size(); ip++){ //loop for all particles
            Qx_2sub[Sub_index[phis_measured_sector[id][ip]]] += cos(ih*phis_measured[id][ip]); 
            Qy_2sub[Sub_index[phis_measured_sector[id][ip]]] += sin(ih*phis_measured[id][ip]); 
          }; // end partilces

          psi_reco_2sub[0] = TMath::ATan2(Qy_2sub[0], Qx_2sub[0] ) / ih;
          psi_reco_2sub[1] = TMath::ATan2(Qy_2sub[1], Qx_2sub[1] ) / ih;
          double dCos_2sub_psi = cos(ih*(psi_reco_2sub[0]-psi_reco_2sub[1]));
          h2sub_Cosndpsi[id][ih]->Fill(dCos_2sub_psi);
        }//end harmonics
      }// 2sub event for V0 done

    } // end detector

    // 2sub evenet doen *****************************************************





    // pTL - EP angle
    for (int id =0; id<ND; id++){
      for (int ih=2; ih<NH; ih++){
        double d_angle = pTL_phi - psi_reco[id][ih];
        if(d_angle< -1*pi) d_angle = d_angle + 2*pi;
        if(d_angle> pi) d_angle = d_angle - 2*pi;    // trick to make hist -pi~pi
        hEP_pTLphi[id][ih]->Fill(d_angle); 
      }
    }

    // calculate each type
    for (int it=0; it<NType; it++){

      int i_Target = Combi_Table[it][kTarget];
      int i_EPSource = Combi_Table[it][kEPSource];

      for( UInt_t ip=0; ip<phis[i_Target].size(); ip++){	
        for( int ih=2; ih<NH; ih++){
          if( i_Target == i_EPSource ){ //remove self-correlation
            psi_reco[i_EPSource][ih]=
                TMath::ATan2(Qy[i_EPSource][ih]-sin(ih*phis_measured[i_EPSource][ip]),
                             Qx[i_EPSource][ih]-cos(ih*phis_measured[i_EPSource][ip]))/ih ;
          }
          double dphi = phis_measured[i_Target][ip] - psi_reco[i_EPSource][ih];
          double real_dphi = phis_measured[i_Target][ip] - psi_in[ih];
          hEPCosndPhi[it][ih]->Fill(cos(ih*dphi));
          hRealEPCosndPhi[it][ih]->Fill(cos(ih*real_dphi));

          hEPVn_Trk[it][ih]->Fill(cos(ih*dphi));
          hRealEPVn_Trk[it][ih]->Fill(cos(ih*real_dphi));

          if(dphi< -1*pi) dphi = dphi + 2*pi;
          if(dphi> pi) dphi = dphi - 2*pi;    // trick to make hist -pi~pi
          if((int)ip <= nPar[it]) hEPdPhi_BulkOnly[it][ih]->Fill(dphi);
          if((int)ip > nPar[it])  hEPdPhi_JetOnly[it][ih]->Fill(dphi);
        }
      }
    }  // calculated done

    delete dNdphi_bulk;


    //   hphi_jet->Write(); hphi_jet->Reset();
    //   hphi_jetpt->Write(); hphi_jetpt->Reset();

    for(int ih=2; ih<NH; ih++){
      for(int it=0; it<NType; it++){
        hEPVn_Evt[it][ih]->Fill( hEPCosndPhi[it][ih]->GetMean() );
        hEPCosndPhi[it][ih]->Reset();

        hRealEPVn_Evt[it][ih]->Fill( hRealEPCosndPhi[it][ih]->GetMean() );
        hRealEPCosndPhi[it][ih]->Reset();
      }
    }
  } // end of the 1 events



  cout <<"Total number of event processed = "<< numberEvents << endl;
  cout <<"Writing histograms into a root file"<<endl;
  fout->cd();
  
  for (int id=0; id<ND; id++){
    hphi[id]->Write();
    hphi_m[id]->Write();
    hJet_phi[id]->Write();
    hDeltaphi_measure[id]->Write();
    hChargedPtDetector[id]->Write();
    hNjetparticles[id]->Write();
  }
  for( int ih=2;ih<NH;ih++ ){
    for( int id=0;id<ND;id++ ){
      hPsi[id][ih]->Write();
      hdpsi[id][ih]->Write();
      hCosndpsi[id][ih]->Write();
      h2sub_Cosndpsi[id][ih]->Write();
//      hEPResolution[id][ih]->Write();
      hJetEPdPhi[id][ih]->Write();
      hJetEPCosdPhi[id][ih]->Write();
      hEP_pTLphi[id][ih]->Write();
//      hEPResolution_2sub[id][ih]->Write();
//      hEPResolution_3sub[id][ih]->Write();
       }
    for( int it=0;it<NType; it++ ){
      hEPVn_Evt[it][ih]->Write();
      hEPVn_Trk[it][ih]->Write();
      hRealEPVn_Evt[it][ih]->Write();
      hRealEPVn_Trk[it][ih]->Write();
      hEPCosndPhi[it][ih]->Write();
      hRealEPCosndPhi[it][ih]->Write();
      hEPdPhi_JetOnly[it][ih]->Write();
      hEPdPhi_BulkOnly[it][ih]->Write();
    }
  }
  hChargedPt->Write();
  hChargedEta->Write();
  hLeadingpPt->Write();
  hLeadingpEta->Write();



  fout->Close();
  timer.Print();
  cout <<"Sucessfully finished."<<endl;
}
